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Abstract 

We consider the assimilation of Lagrangian data into a primitive 
equations circulation model of the ocean at basin scale. The La- 
grangian data are positions of floats drifting at fixed depth. We aim 
at reconstructing the four- dimensional space-time circulation of the 
ocean. This problem is solved using the four-dimensional variational 
technique and the adjoint method. In this problem the control vector 
is chosen as being the initial state of the dynamical system. The ob- 
served variables, namely the positions of the floats, are expressed as 
a function of the control vector via a nonlinear observation operator. 
This method has been implemented and has the ability to reconstruct 
the main patterns of the oceanic circulation. Moreover it is very ro- 
bust with respect to increase of time-sampling period of observations. 
We have run many twin experiments in order to analyze the sensitiv- 
ity of our method to the number of floats, the time-sampling period 
and the vertical drift level. We compare also the performances of the 
Lagrangian method to that of the classical Eulerian one. Finally we 
study the impact of errors on observations. 

1 Introduction 

The world's oceans play a crucial role in governing the earth's weather and 
climate. Lack of data has been a serious problem in oceanography for a 
long time. Since ten years the number of observations has greatly increased, 



with the availability of satellite altimeter data (ie measurements of the free- 
surface height of the ocean) from Geosat, Topex/Poseidon, Jason and other 
satellites. In addition to these remote-sensing data we have in situ data, 
from scientific ships, surface mooring buoys or Lagrangian drifting buoys. 
Among these observations, Lagrangian data, ie positions of drifting floats, 
play a relevant role for many reasons: firstly their horizontal coverage is very 
wide (the whole Atlantic Ocean, for example), secondly they give information 
about currents, temperature and salinity in depth which are complementary 
to surface information given by satellite altimeters. For these reasons many 
national and international programs are organized to deploy drifting floats 
in the world's oceans. The largest program of this type is Argo (2 055 floats 
on the 13th October 2005, 3 000 planned), whose floats provide also temper- 
ature and salinity profiles. 

There are different types of drifting buoys. In the framework of ocean 
basin-scale localized experiments, oceanographers have datasets from acous- 
tic floats. These floats emit acoustic signals which are recorded by moored 
listening stations, and the floats positions are calculated every six hours by 
triangulation. Large datasets are available especially in the Atlantic Ocean 
(SAMBA, ARCANE- Eurofloat, ACCE experiments). On a larger scale Argo 
floats are deployed in order to provide vertical temperature and salinity pro- 
files. Assimilation of Argo thermohaline data has been successfully inves- 
tigated by Forget (PhD Thesis 2004). Argo floats provide also Lagrangian 
information, which are their positions every ten days. Indeed they drift freely 
at a predetermined parking depth (around 1 000 meters), every ten days they 
descent to begin profiles from greater depth (2 000 meters) then they go back 
to the surface and they record temperature and salinity profiles during as- 
cent. On the surface they transmit data to satellite and they are located 
by GPS. Thus many different floats networks and Lagrangian datasets are 
available. 

In parallel, modeling of the ocean system has greatly improved in both quality 
and realism, and there are many Ocean Global Circulation Models (OGCM), 
like for example the OPA PArallelized Ocean model (see Madec et al 1998). 
A crucial issue for oceanographers is then to take the best advantage of dif- 
ferent types of information included in models to one hand and in various 
observations to the other hand. Data Assimilation (DA) covers all theoretical 
and numerical mathematical methods which allow to blend as optimally as 
possible all sources of information (see the review by Ghil et al 1997 and by 
De Mey 1997). There are two main categories of DA methods: variational 
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methods based on optimal control theory (Lions 1968) and statistical ones 
based on optimal statistical estimation (Jazwinski 1970). Adjoint method 
is the prototype of variational methods, introduced in meteorology by Pe- 
nenko and Obraztsov (1976). Its effective implementation in the framework 
of atmospheric Data Assimilation, namely four dimensional variational as- 
similation (4D-Var), has been studied by Le Dimet and Talagrand (1986, see 
also Talagrand and Courtier 1987). Introduction of 4D-Var in oceanography 
is even more recent (see Thacker and Long 1988, Sheinbaum and Anderson 
1990). The prototype of sequential methods is the Kalman filter, introduced 
in oceanography by Ghil (1989) (see also the reviews by Ghil and Malanotte- 
Rizzoli 1991). 

Assimilation of Lagrangian data is in the pipeline. Kamachi and O'Brien 
(1995) used the adjoint method in a Shallow- Water model with upper-layer 
thickness as control vector. More recently Mead (2004) has implemented a 
variational method based on the use of Lagrangian coordinates for Shallow- 
Water equations. Molcard et al (2003) and Ozgokmen et al (2003) imple- 
mented optimal interpolation (which is a simple sequential method) in a 
reduced-gravity quasi-geostrophic model and in a primitive equations model; 
their method is based on conversion of Lagrangian data into velocity informa- 
tion. Ide, Kuznetsov, Jones and Salman (2002, 2003, 2005) used Extended 
and Ensemble Kalman methods to assimilate Lagrangian data into a Shallow- 
Water model; their method is based on an augmented state vector approach 
which does not require the conversion of the positions into velocity data. 
These teams have used simulated data in the twin experiments approach: 
they don't use real Lagrangian data but idealized observations simulated 
from a known "true state" of the ocean. Beside these studies Assenbaum 
and Reverdin (2005) assimilate real data available during the POMME ex- 
periment, including Argo floats data, into a very high resolution model thanks 
to optimal interpolation. 

Previous works on Lagrangian DA were either based on sequential methods 
or on variational ones into very simple models. In this paper we investigate 
variational assimilation of drifters positions into the high resolution primitive 
equations model OPA. 

The aim of variational assimilation methods is to identify the initial state 
of an evolution problem which minimizes a cost function. This cost func- 
tion represents the difference between observations and their corresponding 
model variables. It is minimized using a gradient descent algorithm. The 
gradient is computed by integration of the adjoint model. Thanks to this 
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formulation there is no need to convert Lagrangian data into velocities data: 
we can use directly the position observations, although they are not variables 
of the ocean model, but nonlinear functions of the state variables. Moreover 
this method is a four dimensional one because the temporal dimension of 
the observations, ie their Lagrangian nature, is taken into account. The cost 
function involves a so-called observation operator, which links the state vari- 
ables (here the velocities) and the observed data (here the positions of drifting 
particles). This operator is nonlinear and consequently the cost function is 
not necessarily convex so we used an incremental method (see Courtier et al 
1994) in order to achieve and accelerate the minimization. 
We implement our method using the primitive equations model OPA. Our 
configuration is an idealized wind-driven mid-latitude box model, which is 
representative of the different processes that are going on in the real mid- 
latitude ocean, as shown by Holland (1978). Then we use the twin experi- 
ments approach. As we have said before, Lagrangian observations are simu- 
lated from a known "true state", so that data are perfectly consistent with 
the model and it ensures there is no systematic bias in the observations. It 
is of course unrealistic, but twin experiments are a necessary first-step to 
validate our method. Indeed in this framework we know exactly the system 
true state and so we are able to quantify the efficiency of our method by com- 
paring assimilated and true states. Moreover it was relevant not to use real 
data for many reasons: firstly, Argo floats have not been launched to provide 
Lagrangian information, the feasibility of exploiting their positions is abso- 
lutely not ensured. Secondly, Lagrangian datasets (from Argo to acoustic 
floats) are very diversified in terms of number of floats, time-sampling period 
of observations and drifting depth and we want to investigate the sensitivity 
of our method to these parameters. In order to take into account difficulties 
of real data (such as drift during ascent and descent for Argo floats or acous- 
tic positioning problems) we also study the impact of errors in observations 
on assimilation efficiency. 

The paper is organized as follows: in section [2] we describe the physical model 
and the Lagrangian simulated data. In section [3] we present the assimilation 
method and its implementation. Some numerical results are given and com- 
mented in section HI We conclude in section |U 
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2 Physical model and Lagrangian data 



2.1 The Primitive Equations of the ocean 

The ocean circulation model used in our study is a Primitive Equations (PE) 
model. These equations are derived from Navier Stokes equations (mass con- 
servation and momentum conservation, included the Coriolis force) coupled 
with a state equation for water density and heat equation, under Boussinesq 
and hydrostatic approximations (for more details see Lions et al 1992 and 
Temam and Ziane 2004). 
These equations are written as 

d t u — bAu + (U.V?)u + wd z u — av + d x p = in ft x (0, tf) 
d t v — bAv + (U.V-zjv + wd z v + au + d y p = 
d z p-gT = 

d t T - bAT + (C/.V 2 )T + wd z T + fw = in Q x (0, t f ) { ' 

w(x, y,z) = - Jq d x u(x, y, z') + d y v(x, y, z') dz' in Q x (0,t f ) 

^ U(t = 0) = U , T(t = 0) = T in Q 

where 

- Q = Q 2 x (0, 1) is the circulation basin, where Q 2 is a regular bounded 
open subset of M 2 , x and y are the horizontal variables and z G (0, 1) is the 
vertical one, (0, tf) is the time interval; 

- U — (u, v) is the horizontal velocity, w is the vertical velocity, T the 
temperature and p the pressure; 

- Uo — (uo,vo) and T the initial conditions; 

- (V2O is the horizontal divergence operator and A = d xx + d yy + d zz the 
3-D Laplace operator; 

- a, b, f, g are physical constants. 
The space boundary conditions are 

d z u = t u , d z v = t v , T = on r t 
u = 0, v = 0, T = onffl\r, (2) 
/jLo d x u + d y v dz = in fi 2 

where r = (r u , r„) is the stationary wind-forcing, dVL is the boundary of Q 
and r t = fi 2 x {-2 = 1} is its top boundary. 
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2.2 Model and configuration 

We are using the OPA ocean circulation model developed by LODYC (see 
Madec et al 1998), in its 8.1 version. OPA is a flexible model and can be 
used either in regional or in global ocean configuration. The prognostic vari- 
ables are the three-dimensional velocity field (u,v,w) and the thermohaline 
variables T and S. Discretization is based on finite differences in space and 
time (leap-frog scheme in time). Various physical choices are available to 
describe ocean physics. 

The characteristics of our configuration follows: 

- The domain is Q = (0,1) x (0, L) x (0,H) (longitude, latitude, depth), 
with I = 2800 km, L = 3600 km and H = 5000 m. It extends from -56° to 
-24° West longitude, and from 22.5° to 47.5° North latitude. 

- The horizontal resolution is 20 km, there are 11 vertical levels, so that 
the number of grid points is 180x140x11 = 277200. 

- The time step is 1200 seconds. 

- The model is purely wind-driven. 

This configuration is a classical eddy-resolving double-gyre circulation. As 
shown by Holland (1978) this model is representative of real mid-latitude 
oceans, where circulation is highly nonlinear, non-stationary and where oceanic 
turbulence is very active. Indeed a very active and unstable mid-latitude 
jet develops at the convergence of the subpolar gyre and the subtropical 
gyre. Non-stationary mesoscale eddies form also along the jet. So this model 
shows dynamically different processes such as large-scale gyres, mid-latitude 
jet, mesoscale eddies and also western boundary currents which interact in a 
complex way. Therefore this configuration is a difficult and interesting situ- 
ation to study Lagrangian DA. 

The model is integrated for 25 years until a statistically steady-state is 
reached, which is our "true state" for the twin experiments. Figure [T] shows 
an instantaneous horizontal velocity field at the surface, on the whole hori- 
zontal grid on the left and on a reduced grid on the right. We can see the 
mid-latitude jet and some mesoscale eddies. 

2.3 Lagrangian data 

Lagrangian data are positions of drifting floats. These floats drift between 
Zq — a and z + a where z is given by the user and a is around 25 meters, so 
that we can assume that the floats drift at fixed depth zq, ie in the horizontal 
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Figure 1: Instantaneous horizontal velocity field of the true state. On the 
left the velocity field at the surface on the whole horizontal grid. On the right 
the velocity field at the surface on a reduced grid centered on the mid-latitude 
jet. Dark grey vectors represent large velocities. 

plane z = Zq (Assenbaum 2005, personal communication). We denote by 
£(t) = (£i>£2)(£) the position of one float at time t in the plane z = z . £(t) 
is the solution of the following differential equation: 



where U = (u, v) is the horizontal velocity of the flow and £o the initial po- 
sition of the float. It is important to notice that the mapping U \— > £, which 
links the variables of the model and the Lagrangian observations is nonlinear. 
In the twin experiments approach, observations are simulated by the model. 
The true initial state of the ocean is given and OPA model computes the 
true velocities of the ocean during a ten-day window. We compute on-line 
perfect observations. 

To do that we integrate numerically the equation ^ using a leapfrog scheme. 
This requires the velocity U along the trajectory of the float (ie out of the 
grid). To achieve this we use the following continuous 2D interpolation 




u(t,at),z ) 



(3) 



'interp([7, (x,y)y of the vector field U at the point (x,y): 



x\ = bJ » vi = [y\ > 

u\ = U (zi, t/i), m 2 = ^ {x\ + 1, yi), 

w 3 = U{x l ,y 1 + 1), m 4 = £/(xi + l,y! + 1), 
interp(f7, (x, ?/)) = m + (u 2 - Ui)(x - x x ) + (u 3 - u x ){y - y x ) 
+(ux -u 2 -u 3 + ua)(x - xi)(y - y x ) 

where [-J denotes the floor function, (xi,yi), {x x + {xi,y\ + 1) and 

(xi + + 1) are the grid points which are the nearest neighbors to (x, y). 
This function is piecewise affine with respect to x and y, continuous with 
respect to (x,y), linear with respect to u. Thus it is not differentiable in 
(x, y) everywhere. More precisely it is not differentiable at (x, y) if and only 
if x = x\ or y — y\. It will be a problem to derive the adjoint code, see 



paragraph |3.3| However it is accurate enough to approximate the solution 
of equation ^ and it is very costly to use a differentiable interpolation. 
Indeed such a method (like cubic splines for example) would compute each 
interpolated value from the whole field u (ie the values of u at every horizontal 
grid point) and we would have to inverse a n by n matrix (where n = 25 200 
is the number of horizontal grid points) at every time step and this is not 
workable. 

Let us denote = (^i,fc, ^2,fc) the horizontal position of the float at time 
U the horizontal velocity of the fluid at time Uk the velocity at point 
and h the time step of the ocean model. The algorithm step is schematically 

£fc = £fc-2 + 2/i Uk-i 
U k = interp(C7, £ fe ) 

The dataset is {£at, £ 2 tv, ^3iv---}, where iV is an integer. The duration between 
two data is thus the product of N by the time-step h of the code {h = 1200 
seconds). We call this the time-sampling period. For example if iV = 72 we 
have one data per float and per day and the time-sampling period is thus 
one day. 

In order to simulate real floats we can add errors to the simulated obser- 
vations. Origins of errors are multiple: for acoustic floats, inaccuracy can 
come from acoustic sources (accuracy of their positioning, clocks accuracy, 
bottom topography - acoustic shadow problem, etc.), floats (listening pe- 
riod accuracy, complexity of the trajectory, technical problem - temporary 
"deafness", etc.) or communications quality. For Argo floats errors are due 
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Figure 2: Trajectories of 2 000 floats drifting at level 4 during ten days. On 
the left, trajectories on the whole horizontal grid at level 4. On the right, 
trajectories on a reduced grid around the mid-latitude jet. Very different 
trajectories are observed: short to long, half-circle or straight. 

to drift during ascent and descent and also to drift at the surface between 
ascent /descent and satellite communication. Errors amplitude is around 3 
to 4 kilometers for acoustic floats (T. Reynaud, private communication) and 

2 to 6 kilometers for Argo floats (M. Assenbaum, private communication). 
Figure [2] represents perfect data simulated by the algorithm with 2 000 floats 
for 10 days at level 4 (1000 meters). 

3 Description and implementation of the vari- 
ational assimilation method 

3.1 Description of the assimilation problem 

Without loss of generality we assume that there is only one assimilated data: 
the dataset is the position = (£i(ti), £2(^1)) in the horizontal plane 

z = z of a single float at a single time t\. We use here the notations 
established by Ide et al (1997). We denote by y° = this data. 

Our problem is to minimize the following cost function with respect to the 
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control vector x: 



J(x) = !||£.M(x)-yT + !l|x-x fe ||! 

jT(x) + uj b (x) {V 

where 

- the control vector x = (uo,vq,Tq) is the initial state vector, 

- the B-norm is calculated thanks to the background error covariance 
matrix B _1 : ||^||| = ^B' 1 ^, 

- x fe is another initial state of the ocean, called the background or first 
guess, which is required to be close to the minimum x, 

- M. is the discrete ocean model and .M(x) is the discrete state vector 
(one value per variable, per grid-point and per time step), 

- Q is the discrete nonlinear observation operator, which links the state of 
the fluid (and especially the horizontal velocity U) with the data, QAi{x) = 
£(t) where £(t) is defined by equation (|3]), 

- |.| is the euclidean norm in R 2 , 

- y° is the observations vector. 

Then J° quantifies the misfit between observations and the state of the 
system, J h represents the distance (in terms of the B-norm) between the 
control vector and the background. It is also a regularization term thanks to 
which the inverse problem of finding the minimum x* becomes well-posed. 
The parameter uj represents the relative weight of the regularization term 
with respect to the observation term and it must be chosen carefully. 



3.2 Numerical variational assimilation : incremental 
4D-Var 

Four dimensional variational assimilation (4D-Var, see Le Dimet and Tala- 
grand 1986) is an iterative numerical method which aims to approximate the 
solution x* of discrete assimilation problems with cost function of type Q. 
In 4D- Var a gradient descent algorithm is used to minimize the cost function, 
the gradient being obtained by solving the discrete adjoint equations. It is an 
efficient method but it is very costly when the direct model and the observa- 
tion operator are not linear, for at least two reasons. Firstly every iteration 
of the adjoint method requires one integration of the full non linear direct 
model and one integration of the adjoint of the linearized model. Secondly 
the cost function is not necessarily convex and the minimization process may 
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converge to a local minimum, or it may take considerable time to converge, 
or may not converge at all. 

Incremental 4D-Var (see Courtier et al 1994) avoids, to some extent, both 
of these problems. In this approach, the nonlinear model is approximated 
by a simplified linear model (called tangent linear model) and the nonlinear 
observation operator is linearized around a reference state. The cost function 
becomes quadratical, it has a unique minimum and this minimum is assumed 
to be close to the one of the full non quadratical cost function. In that case 
the minimization process converges quickly. Moreover this approach takes 
into account weak nonlinearities, because the tangent linear model and the 
adjoint model are updated three or four times. 

Remark 1 The approximation of the full nonlinear model by the tangent 
linear model is called the tangent linear hypothesis (TLH). In our highly 
nonlinear configuration, we have to use a ten-day time-window so that the 
TLH is valid (see section^. 



3.3 Implementation in OPAVAR 

The OPAVAR 8.1 package developed by Weaver et al (2003) includes the 
direct non linear model OPA 8.1 developed by LODYC, the tangent linear 
model, the adjoint model and a minimization module. Weaver has imple- 
mented a preconditioning through the B matrix, via the change of variables 
5w = B _1 / 2 5x, following the method introduced by Courtier et al (1994). 
The observation operators of OPAVAR 8.1 are interpolation and projection 
operators. To assimilate Lagrangian data we have implemented the non linear 



observation operator (see section 2.3), its linearization around the reference 
trajectory and the adjoint of the linear observation operator. 
To obtain the tangent (and adjoint) codes of the discrete observation operator 
we use the recipes for (hand-coding) adjoint code construction of Talagrand 
(1991) and Giering and Kaminski (1998). The direct and tangent algorithms 
are schematically: 

• Direct code: 



U k = interp(t/, 

where 'interp' is the interpolation function of U at point £ (see section 
231). 



11 



• Linear tangent code: 

f 6£ k = 6£ k - 2 + 2h6U k -. 1 

\ 8U k = mtexp(8U, £ k ) + 5£ fe .<9( ffi)y )interp([/, £ k ) 

where 'c^^interp' is the derivative of the 'interp' function with respect 
to (x,y). The term S^^unterp is specific to Lagrangian data. It leads 
to a slight difficulty, because the function 'interp' is linear with respect 
of U but it is not derivable in (x, y) at points with integer coordinates. 
Thus we have chosen the values of that derivative at these points, using 
finite centered differences. 



4 Numerical results 

In this section we present the results of our numerical experiments. We begin 
with a brief description of our choices. 

Background and time-window width. In these experiments we have 
assimilated only Lagrangian data and we assume that the true initial tem- 
perature and salinity were known. Background and time-window width are 
related because of the incremental formulation: indeed the full nonlinear 
model is linearized around the background over the whole time- window. When 
the background is too different from the true state or the time-window is too 
wide, approximation errors are large ie the tangent linear hypothesis (TLH) 
is not valid any more. So we compute some correlations to choose both of 
them. If we denote x* the true state and M the tangent linear model, we 
can compute the nonlinear and linear perturbations 8\ and 82 '■ 



S 1 =M(x t ) -M(x b ), $2 = M(x* - 



Then we compute (as a function of time) the spatial correlation between Si 
and 82 according to the formula: 

Cor(Si, S 2 ) = , <W (5) 

y/m-(Si) 2 )m-(W 

where (X) is the spatial mean of X. The closer to 1 the correlation is, 
the better the adequacy between the fields is. Table [T] shows correlation 
at the end of the time-window between linear and nonlinear perturbations 
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Table 1: Background and time-window width choices: spatial correlation 
between nonlinear and linear perturbations at the end of the time-window, 
according to formula The background is the state of the ocean 10 days 
or 1 month before the true initial state. 



time-window width 


background 


correlation 


10 days 


10 days 


0.80 


10 days 


1 month 


0.67 


20 days 


10 days 


0.50 


20 days 


1 month 


0.42 



for different time-window widths (10 or 20 days) and different background 
choices (the state of the ocean 10 days or 1 month before the true initial 
state). We can see that the TLH is not valid with a 20-day window. In 
the sequel we use a 10-day time-window and the state of the ocean ten days 
before the true one as a background state, as in this context the TLH is valid. 
We ran the model with the background as initial state and we obtained a 
"without-assimilation" state, called background in the sequel. It will be 
compared to the assimilated state in order to quantify the efficiency of the 
assimilation process. 

The B matrix. The choice of the B matrix is crucial because of its dual 
purpose (preconditioning and regularization) . Firstly we have tested very 
simple matrices (identity, energy weights) and the results were quite bad: the 
convergence was very slow and the analysis increments were very noisy. These 
matrices are used in OPAVAR only for debugging purpose, with extremely 
idealized assimilation context, for example when the whole state vector is 
observed ie with data everywhere. In order to obtain smoother increments 
and to accelerate the convergence, we used then the diffusion filter method 
(see Weaver and Courtier 2001) which gives good results. 

Diagnostics. Our diagnostics are based on RMS error between the true 
velocity and the assimilated one, compared with the RMS error between the 
true velocity and the background one. The RMS error is plotted as a function 
of time or of the vertical level or of another parameter. For example, we have 
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the following formula for the time-dependent RMS error: 

, .* /Eij,J^J;M) -^,j,M)IV/ 2 ... 

error(M,t = - — — . — - , - — — (6) 

where u t is the true state, u the assimilated state (or the background), t is 
the time and k) a grid point, where are the horizontal coordinates 
and k the vertical one. 

We made the following experiments: first experiment and diagnostics in sec- 



tion 4.1, sensitivity to the floats network parameters (time sampling of the 



position measurements, number of floats, drifting level, coupled impact of 



number and time sampling) in section 4.2, comparison with another varia- 
tional method in section 14.31 

4.1 First experiment 

We present here the results of a typical experiment. There are 3 000 floats 
drifting at level 4 in the ocean for 10 days. The Lagrangian data are collected 
once a day. Thus the total amount of data is 2 x 3 000 x 10 = 60 000. 
Figure [3] (on the left) shows the RMS error of the experiment as function of 
time, according to formula (|6]). We have put the error for the background 
(= without assimilation state) on the same plot. 

Figure [3] (on the right) shows the total RMS error as a function of the vertical 
level (where 1 represents the surface and 10 the bottom), according to the 
formula: 

, j \ &W».i.M)-«(».iiM)IV/ J 

error (u, k) = = , 7 

V 22i,j,t\ u tlh3,k,t)\ 2 J 

We can see that the error with assimilation is twice lower than without. 
Moreover the assimilation process improves every vertical level and not only 
the 4th one. 

Remark 2 We can notice that the RMS error at the beginning are quite 
large. It is explained by the following fact: the relative weight of the reg- 



ularization term J (see section 3.1) has to be large enough to ensure the 



convergence of the minimization process, thus the assimilated state is a com- 
promise between background and observations. 

Figure [1] shows the horizontal velocity field U = (u,v) at level 1 at the final 
time. We can notice that the main patterns as the mid-latitude jet and the 
bigger eddies are quite similar to the true ones. 
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Figure 3: First experiment: u+v RMS errors corresponding to the assimila- 
tion of the positions (sampled 4 times a day) of 3 000 floats drifting at level 
4. On the left, RMS error as a function of time. On the right, RMS error as 
a function of the vertical level. For reference the error without assimilation 
(background) is also displayed. 



true (u,v) 



assimilated (u,v) 




background (u,v) 




Figure 4: First experiment: horizontal surface velocity field at the final 
time. The true field is displayed on the left. The field corresponding to the 
assimilation of floats positions is displayed on the middle. For reference the 
field obtained without assimilation (background) is displayed on the right. 
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Table 2: RMS errors for u+v (in %) at different instants, during a 30-day 
experiment, with and without assimilation of the positions of 1 000 floats 
drifting at level 4, sampled once a day. 



Experiment 


t = 


t = 9 days 


t = 19 days 


t = 29 days 


Without Assimilation 


70.8 


61.8 


62.3 


59.8 


With Assimilation 


55.4 


34.9 


21.7 


13.6 



Longer experiments. We have seen before that the TLH is not valid 
for a 20-day window, so that we cannot use incremental 4D-Var for longer 
windows. However we can restart the assimilation process over the next 10- 
day window in order to do longer experiments: the new background (at the 
beginning of the new window) is the previous assimilated state (at the end of 
the previous window), so that the new background carries information from 
the previous assimilation process. Table [2] shows the relative RMS errors 
([7]) for u + v of a 30-day experiment at different times. We can see that 
error with assimilation is lower than 15% at the end of the window, ie it is 
less than one fourth of the error without assimilation. This is a very good 
result: 3 succesive assimilation processes enable to reconstruct a very good 
approximation of the true state. 



4.2 Sensitivity to the floats network parameters 

In operational oceanography, sensitivity analysis is central to the observa- 
tional network design. Indeed, in situ and remote observation instruments 
are very expensive (e.g. one Argo floats costs 15 000$) and they must be op- 
timally used. Ngodock (PhD thesis 1996) shows that second order analysis 
(ie the derivation of the optimality system or second order adjoint system, 
see also Wang et al 1992) enables to analyze the sensitivity of the 4D-Var 
assimilation system to the design of the observational network. Second order 
adjoint information (see also the review paper by Le Dimet et al 2002) is 
actually central to adaptive observation network and observation targeting 
issues. However it requires the storage of model, tangent and adjoint trajec- 
tories, so that it is not workable in OPAVAR at the present time because of 
computer memory limitations. 

So we have performed a lot of experiments to analyze the sensitivity to various 
parameters of our assimilation process. Indeed in the ocean the network's 
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parameters can widely change, from Argo floats (1000-2000 meters depth, 
one data per 10 days) to acoustic floats (various depth, time sampling period 
around 6 hours) or drifters in the upper ocean (near surface, time sampling 
period very short)... 

Here are the parameters that we consider: 

- the time sampling period, varying from 6 hours to 10 days, 

- the number of floats, varying from 300 to 3000, 

- the vertical level of drift, varying from 1 (surface) to 10 (bottom), 

We analyze also the coupled effect of the number and the time sampling 
period. 



4.2.1 Sensitivity to the time sampling period. 

The framework of this experiment is the following: we performed seven differ- 
ent experiments with exactly the same initial conditions, namely 3 000 floats 
at level 4. The only difference in these experiments is the time sampling pe- 
riod, which will be denoted shortly by TSP in the sequel. The experiments 
are denoted by TSP-xxx where xxx is the time sampling period, in hours 
(6 or 12h) or in days (1, 2, 3, 5 or lOd). Figure [5] shows the RMS error as 
a function of time for each TSP experiment, except (for readability) exper- 
iments TSP-6h, TSP-12h and TSP-2d. It shows also the total RMS error 
as a function of the time sampling period. We can see that our method is 
robust with respect to the increase of the time sampling period. This is very 
encouraging. Indeed every prior study is very sensitive to the TSP and shows 
quite bad results when the TSP is larger than 2 or 3 days (see Molcard et al 



2003, Mead 2005 and section 4.3). Our method does not show this sensitiv- 
ity, velocities are quite well reconstructed even when the TSP is large and 
especially with a 10-day period, which is a very positive result. 



4.2.2 Sensitivity to the number of floats. 

We perform five experiments with varying number of floats drifting at the 
same vertical level (4) and with positions sampled with the same period (6 
hours). The experiments are denoted by NUM-xxx, where xxx is the number 
of floats. Figure [6] shows the RMS error as a function of time and the total 
RMS error as a function of the number of floats for each experiment. We 
can see that the number of floats has great influence on the results. Under a 
minimal number (1 000) the velocities are badly reconstructed, undoubtedly 
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Figure 5: Sensitivity to the time sampling period: u+v RMS errors corre- 
sponding to assimilation of 3 000 floats' positions with different time-sampling 
periods of observation. On the left: error as a function of time for each 
'TSP' experiment and for the background (without assimilation reference 
state). On the right: final error as a function of the TSP; the error without 
assimilation is also displayed. 
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Figure 6: Sensitivity to the number of floats: u+v RMS errors corresponding 
to assimilation of 300 to 3 000 floats' positions. On the left: error as a 
function of time for each 'NUM' experiment and for the background (without 
assimilation). On the right: final error as a function of the number of floats; 
the error without assimilation is also displayed. 

because there is not enough information to constrain the flow. However, when 
we perform longer experiments, we get satisfactory results for small numbers 
like 500 and 300. For a ten-day window the results are optimal with a 1 000 
floats network and they don't improve with higher numbers. Obviously the 
information becomes redundant and it is useless to add floats. The associated 
density is one float per 10 000 km 2 , which is ten times more than the planned 
Argo density (namely around 100 floats in our configuration). Even if we 
perform longer experiments, the Argo density is too small to constrain the 
velocity field. It is more appropriate to use Lagrangian data from localized 
experiments such as acoustic floats launchings in the Atlantic Ocean (like 
e.g. SAMBA), whose floats densities are higher. 

4.2.3 Sensitivity to the vertical drift level. 

Again we perform seven experiments with 3 000 floats drifting at varying 
vertical level and fixed TSP (6 hours). As usually we denote by LEV-x the 
experiment involving floats at level x. Figure [7] shows the RMS error as a 
function of time for upper levels on the left and lower levels on the right. 
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Figure 7: Sensitivity to the vertical drift level: time-evolution of the u+v 
RMS errors corresponding to the assimilation of positions of floats drifting at 
different depths. On the left: error for upper levels experiments (above 1 000 
meters) and for the background (without assimilation). On the right: error 
for lower levels experiments (below 1 000 meters) and for the background 
(without assimilation). 

Again the results are very sensitive to the position of the floats. The three 
best levels are 3, 4 and 5, ie the intermediate levels. From a physical point of 
view it is coherent because the information propagates vertically with a finite 
velocity so that very upper (1, 2) and very lower (7 to 10) levels are penalized. 
Moreover upper levels (1 to 4) are the most energetic ones (from the kinetic 
turbulent energy point of view), quasi ten times more than the lower ones 
(levels 5 to 10), it seems quite natural that the best results are obtained with 
floats drifting at level 4 which is both intermediate and energetic. 

4.2.4 Coupled impact of number of floats and time sampling pe- 
riod. 

Here we look at the coupled effect of varying number of floats and vary- 
ing TSP, for example in order to answer the following question: is the total 
number of data an important variable to measure the efficiency of the assim- 
ilation? So we perform nine experiments, denoted by nnn-xxx where nnn is 
the number of floats and xxx is the TSP. These experiments and their final 
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Table 3: Coupled impact of number of floats and TSP: Final RMS error cor- 
responding to assimilation experiments with 500 to 2 000 floats and positions 
sampled every 1 to 5 days. Total number of observations is also given. The 
"background experiment" results are also shown for reference. 



Experiment 


Total number of data 


Final Error (%) 


500-5D 


1000 


46.6 


500-3D 


1500 


44.0 


500-1D 


5 000 


44.0 


2000-5D 


4000 


37.1 


2000-3D 


6000 


36.8 


2000-1D 


20000 


35.9 


1000-5D 


2 000 


35.1 


1000-3D 


3000 


35.1 


1000-1D 


10 000 


34.8 


Background 


no data 


60.8 



RMS error are described in Table [3] Figure [8] represent the RMS error as 
a function of time for the 500-xxx experiments on the left, 1000-xxx in the 
middle and 2000-xxx on the right with the same scale on the axis of ordi- 
nates. The results are complementary to the precedent experiments. Indeed 
we can see that 1 000 is an optimal number for this configuration whatever 
the TSP and that our method is stable with respect to large TSP whatever 
the number of floats. Thus we can conclude that, in our configuration, it 
seems optimal to launch around 1 000 floats and that the TSP can be chosen 
quite large. 



4.3 Comparison with the "Eulerian" method 

A classical method in oceanography is to assimilate the velocity observations 
deduced from the Lagrangian data according to the following finite differences 
formula: 

£l(4+l) -6(4) ^ (£ (+ \ £ t+ \ 4- \ 

£2(*fc+l) -6(4) ( C ( 4- \ £ (4- \ 4 \ 

— 7 7 ~ ^(£i(4),6(4),2o,4) 

£fc+i — £fc 

Then the velocity data are treated as Eulerian data (measured at non-fixed 
points). We implement this method in the 4D-Var framework. The obser- 
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Figure 8: Coupled impact of number and time sampling period: time- 
evolution of the u+v RMS errors corresponding to assimilation experiments 
with 500 to 2 000 floats and positions sampled every 1 to 5 days: experiment 
with 500 floats on the left, 1000 in the middle and 2 000 on the right. For 
reference, the background error is also displayed on each plot. 



vation operator is much easier to write (and to differentiate and transpose) 
because it is an interpolation at the points of the true (fixed) floats trajec- 
tories. We compare the results for this method said "Eulerian" and for our 
"Lagrangian" one. Experiments have the same characteristics (3 000 floats 
and varying TSP), their names are LAG-xxx or EUL-xxx with xxx the TSP, 
where xxx is the TSP. 

Figures [9] and 10 represent the RMS error as a function of time and vertical 
level. We can see that the "Eulerian" approach is slightly better than the 
"Lagrangian" one when the TSP is small (one day), moreover its computa- 
tion time is 10% lower. Indeed the TSP is small enough so that the formula 
Q is a very good approximation: the displacement vector between two suc- 
cessive positions is quasi tangent to the trajectory (ie quasi collinear with 
the velocity vector). However we can see that the error for the "Lagrangian" 
method is more homogeneous as a function of the vertical level. For larger 
TSP (3 days or more) the "Lagrangian" method is obviously better than the 
"Eulerian" one: the approximation formula ^ is not valid any more. We 
can see that our method is able to extract information from the positions 
data even if the TSP is large. 
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Figure 9: Comparison Eulerian/Lagrangian: time-evolution of the RMS error 
for u+v corresponding to the assimilation of 3 000 floats with different TSP. 
On the left, errors for the Lagrangian and Eulerian methods with small TSP. 
On the right, errors for the Lagrangian and Eulerian methods with large 
TSP. For reference, the background error is also displayed on each plot. 



4.4 Assimilation of noisy observations 

In order to deal with real data issues, a necessary first-step is to study the 
impact of observation errors in the twin experiments framework. To do that, 
we simulate as previously perfect data from the "true state" with 1 000 floats 
drifting at level 4, their positions being sampled once a day. Then we add a 
random Gaussian noise to the computed positions. In the sequel, the word 
"error" represents the amplitude of the noise. As we said before, real errors 
are about 2 to 6 kilometers. However our system is idealized so we study 
the impact of errors up to 20 kilometers. The total displacement of one float 
between initial and final positions is around 25 kilometers (in steady regions) 
to 90 km (in the mid-latitude jet region), so that a 10 to 20-km noise is sig- 
nificant for most of the floats. 



Figure 11 represents RMS errors as a function of time (on the left) for exper- 
iments with to 20km errors and for the background (without assimilation). 
On the right we plot RMS errors as a function of observation error amplitude. 
The RMS error is very stable with increasing noise amplitude: our method is 
able to extract information even when the error amplitude is not negligible 
with respect to floats displacement. 
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Figure 10: Comparison Eulerian/Lagrangian: final RMS error for u+v as 
a function of the vertical level, corresponding to the assimilation of 3 000 
floats with different TSP. On the left, errors for the Lagrangian and Eulerian 
methods with small TSP. On the right, errors for the Lagrangian and Eulerian 
methods with large TSP. For reference, the background error is also displayed 
on each plot. 
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Figure 11: Impact of observation errors: u+v RMS errors corresponding to 
assimilation of noisy observations. On the left: error as a function of time for 
experiments with noise amplitude from 1km to 20km. For reference, results 
for experiments without noise (0km) and without assimilation (background) 
are also displayed. On the right: final error as a function of the amplitude 
of the noise; the error without assimilation is also displayed. 



Figure 12 shows the evolution of the cost function value and its gradient 
during the assimilation process. The abscissa represents the total number of 
iterations. Here we perform four outer loops: in each outer loop we perform 
ten inner minimization loops, in which the cost function is minimized, as we 
can see on the left. We can see that the gradient norm decreases and the 
cost function converges even in the presence of noise in data. 



5 Conclusion 

This paper shows that the problem of assimilating Lagrangian data can be 
solved by a variational adjoint method into a realistic primitive equations 
ocean model. We have implemented a Lagrangian method which takes into 
account the four dimensional (space and time) nature of the observations: 
RMS errors with assimilation are twice lower than without and the main 
patterns of the fluid flow are well identified at each vertical level, although 
the floats drift at a single determined level. 
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Figure 12: Evolution of the cost function and its gradient's norm during the 
assimilation of noisy data. On the left, evolution of the cost function for 
experiments with noise amplitude from 1km to 20km. For reference, the cost 
for experiment without noise (0km) is also plotted. On the right, evolution 
of the gradient norm, shown on a logarithmic scale, for experiments with and 
without noise in data. 

We have tested the sensitivity of our method to the characteristics of the 
dataset. It is very sensitive to the vertical drift level, and the best results 
are obtained for intermediates ones, especially level 4 (around 1000 meters 
depth). It is also very sensitive to the number of floats, but the more is not 
the better, it seems useless to launch more than 1 000 floats in our config- 
uration. It is very robust with respect to the increase of the time-sampling 
period, up to ten days. 

We have compared our Lagrangian method to the Eulerian one, which con- 
sists in interpreting Lagrangian data as velocity information. When the 
time-sampling period of the observations is one day or less, the Eulerian 
method performs slightly better, but the transfer of information to lower 
levels is better achieved by the Lagrangian one. When this period is larger 
than two or three days, the Lagrangian method performs much better than 
the Eulerian one. 

We also studied the impact of errors on observation: the reconstruction of 

the velocities is well achieved even with a large noise in data. 

Also the performances of this method have been assessed in the framework 



26 



of the twin experiments approach. The next step would be to use real data 
and to deal with problems such as trajectories modelling and model error. 
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